 function w = newton(x,y,z)

    w = [];
    dv = dif_divididas(x,y);

    for i = 1:length(z)
        w(i) = pol_newton(x, dv, z(i));
    endfor

endfunction

function d = dif_divididas(puntosX, puntosY)
    N = length(puntosY);
    for i = [1:N];
        d(i) = puntosY(i);
    endfor
    
    for j = [2:N]
		for i = [N:-1:j]
            d(i) = (d(i)-d(i-1)) / (puntosX(i)-puntosX(i-j+1));
        endfor
    endfor        
endfunction

function v = pol_newton(puntosX, dv, x)
    N = length(puntosX);
    s = 0;

    for i = [1:N]
        j = 1;
        p = 1;
        while (j <= i-1)
            p = p * (x-puntosX(j));
            j = j+1;
        endwhile
        s = s + dv(i)*p;
    endfor
    v = s;
endfunction

